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We predict the structural interaction of crystalline solid-melt interfaces using amplitude equations 
which are derived from classical density functional theory or phase-field-crystal modeling. The solid 
ordering decays exponentially on the scale of the interface thickness at solid-melt interfaces; the 
overlap of two such profiles leads to a short range interaction, which is mainly carried by the 
longest-range density waves, which can facilitate grain boundary premelting. We calculate the tail 
of these interactions, depending on the relative translation of the two crystals fully analytically and 
predict the interaction potential, and compare it to numerical simulations. For grain boundaries 
the interaction is predicted to decay twice faster as for two crystals without misorientation. 
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I. INTRODUCTION 

Grain boundaries (GBs) and interfaces in general have 
a strong influence on mechanical behavior and other ma- 
terials properties. Therefore, they have been widely stud- 
ied both experimentally^ and computationall}Elfor a long 
time. At high temperatures close to the melting point, 
GBs can display pronounced disorder, even leading to 
the formation of nanometer-scale intergranular films with 
liquid-like properties. The formation of those films below 
the bulk melting point, known as GB premelting, lead to 
catastrophic materials failure, initiated by an enormous 
reduction of the shear resistance. This phenomenon is 
of interest for predicting the formation of solidification 
defects associated with the formation of those intergran- 
ular films, which can lead to hot cracking during the late 
stages of solidification^ ^, and more generally for under- 
standing the microstructure and mechanical behavior of 
structural alloys at high homologous temperature. 

Over the years, there have been many experiment aP^^^ 
and theoretical investigations of GB premelting. The lat- 
ter include discrete lattice models^^ and molecular dy- 
namics (MD) simulationJI^, as well as conventional 
phase-field modelj^^^ HI^ whi ch either exploit an ori ent a- 
tional order paramete j^^ * ^^ * or multiple phase-field^^^^ 
to distinguish between grains, and the phase-field-crystal 
(PFC) method^^^^, which resolves the crystal density 
field on an atomic scale and hence naturally models crys- 
tal defects such as isolated dislocations and GBs. 

The determination of the premelted layer width W and 
the quantification of the fundamental forces that control 
this width are of striking interest for GB premelting. Ex- 
perimentally, these issues are difficult to address. Obser- 
vations to date support the existence of a nanometer- 
thick premelted layer in pure materials a few degrees be- 
low the bulk melting point and there is more ample evi- 
dence for premelting in alloys. Modeling activities based 
on PFG^ and MD simulation^— allow characteriza- 
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FIG. 1. Sketch of the atomic density profile at two adjacent 
solid-melt interfaces and the corresponding amplitude. 



tion of the structural forces underlying GB premelting. 

These forces are quantified by the introduction of the 
"disjoining potential" T^(VF), which is defined via the ex- 
cess Gibbs free-energy per unit of grain boundary area 

Gexc(W, T) = AG{T)W + 2jsi + V{W), (1) 

where AG = Gs — Gi is the bulk Gibbs free-energy dif- 
ference between liquid {Gi{T)) and solid (Gs(T)) and 
is the solid-liquid interfacial free-energy. Based on this 
definition, represents the part of this excess due to 

the overlap of crystal density waves from the two grains 
on each side of the GB (see Fig. [T]). The ordering of the 
solid phases extends also into the melt on the range of the 
interface thickness ^, thus the crystals start to interact 
with each other as soon as their separation W is of the 
order of the interface thickness. Depending on the align- 
ment of the crystals their structures may match — which 
leads to an attractive interaction — or are locally dis- 
placed such that the energy of the system is increased by 
the overlap, which leads to a repulsive interaction. Hence 
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the derivative —dV{W)/dW expresses the force between 
crystal- melt interfaces due to this overlap, which can be 
either repulsive or attractive depending on whether the 
sign of —dV{W)/dW is positive or negative, respectively. 
So far, there is little analytical knowledge on the short 
range contributions to these forces, with the exception 
of phase-field models^^, which are based on phenomeno- 
logical models that do not consider atomic structures, 
dislocation formation and elastic interactions. The pur- 
pose of the present article is therefore to gain analytical 
insights into the nature of these forces, based on a com- 
plex Ginzburg-Landau description. Additionally, 1/(1^) 
also contains attractive contributions due to London dis- 
persion forces that are neither accounted for in PFC and 
MD simulations, nor in the short range forces analyzed 
here, but play an iniportant role in other systems such 
as ceramic materials^. However, in metallic systems, 
dispersion forces can be estimated to only contribute an 
attractive tail to ^(W) whose magnitude is less than a 
mJ/m^ for W on the nanometer scale. In contrast, MD 
computations of V{W) in pure Ni-^^-^ show that V{W) 
has a magnitude of tens of mJ/m^ for W in this same 
range. 

A major outcome of this work is that the structural 
interaction can be calculated analytically in the case of 
zero misorientation between the grains, which only have a 
translational misfit. The range of the interaction can still 
be computed also for the more general case of misoriented 
crystals. The results therefore offer new insights into 
the phenomenon of GB premelting, as they show which 
quantities and ingredients are essential for the structural 
interactions. This paper therefore complements the nu- 
merical results in Ref. 31, which compare PFC and am- 
plitude equations results with MD data. 

The structure of the article is as follows: In Section 
[IT| the underlying model is summarized, which is then 
used in Section |lll] to investigate the properties of single 
solid-melt interfaces. There, in particular, the decay of 
the density waves into the melt is analyzed, since this 
turns out to be the key parameter for the interaction of 



two solid- melt interfaces. First, in Section |IV| the special 
case of crystals without misorientation but with transla- 
tional misfit is considered, as here the asymptotics of the 
interaction can be calculated fully analytically, which is 
demonstrated for several interface orientations. In Sec- 
tion [V| grain boundaries are considered; although a full 
analytical treatment is not possible here, still the range 
of the interaction can be predicted. 



II. AMPLITUDE EQUATIONS 

From the classical density functional theory of freezing 
(DFT) a functional can be derived which expresses an 
emerging solid phase as density fiuctuations appearing 
from the liquid state, whereas the (time- averaged) den- 
sity is spatially constant in the melt phase^^^*^. To this 
end the spatial variations of the density field ^'^(r), are 



expanded as a sum of density waves 

N 



(2) 



where each ki is one of the N different reciprocal lat- 
tice vectors (RLVs) and u^^"^ are their associated am- 
plitudes. In the liquid phase, where the time averaged 
density is spatially constant, the amplitudes vanish, and 
in an undistorted solid phase they all attain the same 
constant value, u^^"^ = Ug. The associated free energy 
deviation from the liquid state is 
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where the function /({i^^-^^}, T) contains the higher order 
nonlinear terms in the amplitudes u'^^^ and an explicit de- 
pendence on the temperature T. Furthermore, no is the 
density in the liquid state and C{r) is the direct correla- 
tion function with Fourier transform 



C{q) = no / dfC{r)ex.-p{—ik-f) 



with r = \r\ and q 
ture factor by 



(4) 



\k\; it is related to the liquid struc- 
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Here, all quantities are evaluated at the (first) peak of the 
structure factor go- Unlike DFT, where a large number of 
modes is required to obtain sharp peaks around atomic 
positions, the simpler free energy allows for a truncation 
of this sum to a small set of reciprocal lattice vectors. 
Various methods have been developed to change the ker- 
nel of the free energy in order to stabilize a variety of two 
and three dimensional periodic and crystal structures. 
Here we focus on bcc structures, therefore we restrict the 
summation to the N = 12 principal reciprocal vectors 



[110], [101], [Oil], [110], [101], [Oil] 
[110], [101], [Oil], [110], [101], [Oil]. 



(6) 



Notice that by the condition of having a real density field 
i/j the N complex amplitudes are not independent but 
are complex conjugate (denoted by a star) for antipar- 
allel RLVs. Therefore, we restrict the description to the 
first row of these RLVs and use only independent 
complex fields u^^\ 

The differential operator is given b}^^^^^ 



2go 



(7) 



where the vectors k^^"^ are the normalized principal RLVs. 
The second term in the operator preserves the rotational 
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invariance of the equations and is related to the use of 
the nonhnear strain tensor in elasticity. Below we refer 
to this second term as the higher order correction term in 
the box operator, since it is vanishingly small for rough 
interfaces. 

For most of the present analysis the precise form of 
the higher order nonlinearities is not important, as they 
enter the expressions for the interface interaction only 
as pref actors in terms of matching constants. However, 
to complete the model, we use here amplitude equations 
which are derived via a multiscale expansion from the 
three-dimensional phase-field-crystal modet^. 



2 ^ 
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Here, e is used as a small parameter in the regime of the 
phase diagram which describes the coexistence between 
the bcc and the homogeneous (melt) phase. The param- 
eter e is related to the physical parameters via 

with another (small) dimensionless parameter e 

24 

6 = 



(10) 



which turns out to be more useful in the context of the 
amplitude equations; it characterizes the ratio between 
the square of the atomic spacing and the solid-liquid in- 
terface thickness. 

In equilibrium the chemical potential 



(11) 



is spatially constant. A detailed derivation of the ampli- 
tude equations, which describe the evolution of the fields 
u^^^ has been given in Refs. 38 and 40, and therefore we 
only give the resulting expressions here. The evolution 
equations can be derived from a free energy functional, 
which reads 

■ Ar/2 



1 

12 



AT/ 2 




AT/ 2 
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+2Aiio^oii^oii^iio + 2A*jqAq;li^oii^iio 



+274oiiA^Qj74ioiAoii + 274q-|^j74xoi^ioi^ 



oil 



~g (^011^101^110 +^011^101^110 +^011^110^101 
+^011^110^101 + ^011^110^101 + ^011^110^101 
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Here, we have introduced rescaled amplitudes 

^(j) =^^0-)/^^^. (13) 

Also, we have introduced a dimensionless length scale 

R = e^/^^or. (14) 
In these rescaled coordinates the box operator becomes 

-1/2 ^ 



and the prefactor Fq is 



C'\qo)qo'ul~e 
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(15) 
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Finally, the thermal tilt Ft is added phenomenologically, 

■ Tm 



dRL^, 



-0({^(^)}) (17) 



to favor either the solid or liquid state. Here, Tm is the 
melting temperature, L the latent heat and (j) is an "order 
parameter" which discriminates between solid and liquid. 



N/2 



with 



h{(f)) = (/)2(3 - 2<P). 



(18) 



(19) 



Alternatively, the tilt can be chosen such that it repro- 
duces the PFC results, and then the coupling function is 
chosen to be 



N/2 



(20) 



We note that this expression is invariant under elastic 
deformations and lattice rotations, which affect the com- 
plex phases of the amplitudes. In the original DFT for- 
mulation the coupling therefore reads 



Ft = L- 



T-Ti 



M 



Tm 



i=i 



(21) 



We note that the two above coupling functions (18) 



and (20) are substantially different: The first is quar- 



tic in the amplitudes variations in the solid and liquid 
state, and therefore the minima of the functional remain 
at A^^ = and A^^ = 1 for T 7^ Tm- This is not the 
case for the second coupling, which is linear in the am- 
plitudes; therefore here the bulk values depend on the 
temperature. We will discuss the implications of these 
two different couplings in Appendix [Aj 
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Thermodynamic equilibrium corresponds to a station- 
ary state of the free energy functional, and we use relax- 
ation dynamics according to 



dt 



5F 



(22) 



Since we focus here exclusively on static situations, the 
choice of the kinetic coefficients Kj is arbitrary. 

This description predicts the correct anisotropy of sur- 
faces energies^^ and elastic properties and contains nat- 
urally the linear theory of elasticit}' . The form of these 
nonlinearities depends slightly on the underlying model: 
Above it is given for a PFC model, and there are some 
differences if these terms are derived from DFT using an 
equal weight ansatz. However, the differences are small 
and lead e.g. only to tiny changes of the anisotropy of 
the surface energy, as had been investigated in Ref. [391 
Also, we point out that — as will be shown in the fol- 
lowing sections — the higher order nonlinearities do not 
contribute to the short-range interaction tail for shifted 
crystals. 

Finally, we note that this amplitude equations model 
for crystals is conceptually close to theories of supercon- 
ductivity and pattern formation in hydrodynamics^. 



III. SOLID-LIQUID INTERFACES 

Properties of solid- melt interfaces, in particular inter- 
facial energies and their anisotropy, were discussed in de- 
tail in Ref. |38] and |39| Here we concentrate on specific 
properties that are relevant for the understanding of in- 
terface interactions in the next section. 

In the melt region sufficiently far away from the inter- 
face, the amplitudes of the density waves have decayed 
and can be well described by the linearized equilibrium 
conditions or, equivalent ly, the free-energy density with 
the local terms up to second order in the amplitudes. It 
is worthwhile to mention that the free energy functionals, 
as derived from PFC and DFT, agree up to this order, 
and therefore the following results are generic. 

Let the interface normal n be the z direction of a 
straight solid-liquid interface, and all density wave ampli- 
tudes depend then only on this coordinate. In the liquid, 
where the solid ordering has decayed almost completely, 
the equilibrium conditions decouple and are given by 



S{qo) 



(23) 



where we ignore for the moment the higher order cor- 
rections of the box operator. For reasons that will be- 
come more obvious later, we denote here derivatives with 
respect to 2; by a dot. Although we consider a three- 
dimensional situation, the amplitudes depend only on 
z. We focus here on stationary states, therefore time- 
derivatives do not appear. Obviously, (stationary) coex- 
istence between solid and liquid bulk phases with a planar 
interface is only possible for T = Tm- 



The general solution of this linearized equilibrium con- 
dition is a superposition of two exponentials, 

u''^^ = Cj^in ex.p{-Xjz) + Cj^out ex.p{Xjz) (24) 

with the inverse decay length 

-2 
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We also define the characteristic scale Aq ~ 1/^ 



Ao 



1/2 



(25) 



(26) 



Since we consider only a single interface, with the solid 
phase being located at z ^ — 00, the growing exponential 
cannot be present {cj^out = 0). We note that a shift of 
the interface by Az leads to a change of the remaining 
exponential prefactor Cj^in by a factor exp(— AjAz). 

Since the problem is one-dimensional, it is straightfor- 
ward and fast to solve the full set of amplitude equations 
(not only in the linearized region) using a real space im- 
plementation via the relaxation scheme (22) at T = Tm, 



until an equilibrium solid-liquid interface is established. 
The grid spacing is chosen to be much smaller than the 
interface thickness, Xjdz <C 1, to obtain results which 
do not depend on the discretization. Corresponding to 
the analytical investigation we do not take into account 
the higher order term in the box operator, thus the equi- 
librium profile is described by a second order ordinary 
differential equation. From the equilibrium profile the 
solid-liquid interfacial energy jsi is also computed. To 
obtain a numerical value for the prefactor Cj^in we have 
to match it to the full solution of the nonlinear problem 
5F/5u^^ = 0. Therefore, we set the origin z = exactly 
at the location of the interface. Since the interface is 
smooth, the position of it requires a precise definition. 
The choice of this measure is not critical, since another 
definition only leads to slightly different values for the 
exponential prefactors, and later on in the following sec- 
tions, to a horizontal shift of the disjoining potential. We 
use an integral measure for the amount of liquid per unit 
area of the interface 
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dzY,[l-h{\u(^^\yu% (27) 



where we use the "coupling function" ( |19[ ), h{x) = 
x'^{3 — 2x) to interpolate between solid and liquid. Notice 
that in the liquid the amplitudes have decayed to zero, 
whereas in the solid all of them have the value 

For a single solid-liquid interface, the amount of liquid 
depends of course on the system size, i.e. the length of the 
integration interval L^. For Lz ^ zq {zq is the interface 
position), W becomes a linear function of L^, W Lz — 
Zq. We can extrapolate this linear function to the value 
0, which then defines the location of the interface, and 
this is shown in Fig. [2] 
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FIG. 2. Determination of the interface position. The sohd 
phase is located in the right, the hquid in the left half of the 
system. The amount of melt is defined through the measure 
(27). For sufficiently large systems, this expression becomes 
asymptotically equal to — , where Lz is the system length 
and zq the interface position. Results are shown here for a 
(100) interface at T = Tm- In the present case, the interface 
is located at zqXq = 21.7. 
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FIG. 3. Matching of the exponential decay of the amplitudes. 
For a (100) surface, the amplitudes group into two classes. 
The matching constants are determined such that the curves 
for the slowest decaying density wave amplitudes (solid curve) 
agree with the exponential solution of the linearized equations 
(dotted line). 



In the next step, we plot the amplitudes as a func- 
tion of the distance from the interface, z — zq. For 
z ^ zq they decay exponentially on the scale 1/Aj, and 
we can determine the exponential prefactors as shown 
in Fig. [3j We note that for a single interface all ampli- 
tudes can be chosen to be purely real (as long as the 
correction from the box operator is suppressed). For a 
(100) surface, only the density waves [110], [110], [101], 
[101] (+ complex conjugates) follow the exponential de- 
cay given above, since for the others the interface normal 
is perpendicular to the principal reciprocal lattice vec- 
tors, /c • n = 0. This means that they decay faster, as 
they are "slaved" by the other fields (see Fig. [3|. As 
explained in Ref. the density waves can therefore be 
grouped into two classes. The matching constant for the 



FIG. 4. Matching of the exponential decay of the ampli- 
tudes. For a (110) surface, the amplitudes group into three 
classes: [110] has the longest range, the second group is 
[101], [101], [Oil], [Oil]. The shortest range density wave is 
related to the wave vector [110]. 



slowly decaying fields is then determined numerically as 

^110^ — ^iio^^ ~ ^101^^ — ^loi^^ ~ 0.165 for the defini- 
tion of the interface position given above (the subscript 
denotes the density wave, the superscript the interface 
normal) . 

The same procedure can be repeated for any other in- 
terface, and in general all matching constant are differ- 
ent (they are only pairwise equal for complex conjugate 
fields). The corresponding plot for a (110) surface is 
shown in Fig. l4j The obtained matching constants are 

^110 — u.iio ana c^^oi — ^loi ~ ^Oll ~ ^oii ~ 
0.372; the remaining field Uiiq is slaved by the others 
and decays faster. 

Finally, for a (310) interface (as a representative for 
an arbitrary interface normal direction), the numerical 

0.18, 



matching gives: Cy^^^ = 0.128, c^q^^^ 

^(310) 
-^110 



.(310) 
-^101 



= 0.425. We note that the fields uqh and 



which are expected to have a decay according to Eq. ([25| , 
turn out to behave differently; they decay more slowly 
than anticipated, because e.g. for uqh a forcing term of 
the structure ~ i^iio'^ioi ^^^^ the cubic terms in the free 
energy functional leads to a longer range of the density 
wave than the anticipated quadratic term; in fact, for the 

given inclination Aq^^^^ > A^^q^^ + ^lol^^ ^^^^ right 
hand side being the inverse decay length of the slaved 
field. In general, it means that also density wave ampli- 
tudes with k ■ h can be slaved by other terms, and 
their decay is then determined by the higher order non- 
linearities. The range of these slaved fields is very short, 
and therefore they do not contribute significantly to the 
interaction potentials derived below. The decay of all 
amplitudes, together with the exponential fits, is shown 
in Fig. [5] 

Let us briefly discuss the relevance of the box opera- 
tor corrections to the preceding results. So far, a single 
straight interface is in principle described by a real den- 
sity wave amplitude, or at least by a constant complex 
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FIG. 5. Matching of the exponential decay of the ampUtudes 
for a (310) sohd-hquid interface. The range of the fastest 
decaying field uqh is not determined by quadratic but cubic 
interactions. 



phase. The box operator explicitly introduces an imagi- 
nary factor, and therefore the amplitudes pick up a small 
and slow oscillatory contribution. Also, the results with- 
out the higher order corrections (as given on the slow 
scale), do not yet depend on the value of e, which is only 
re-introduced when the density profile is reconstructed 
from the amplitudes. With the correction terms, e ap- 
pears also explicitly in the amplitude equations. The 
linearized equations in the liquid region become 



S{qo) 



+ 



Jc"(9o)n?««) = o, 



(28) 



and in one dimension 



' ^ ' 2|fcO)| 



U) 



(29) 



The general solution of the linearized problem then be- 
comes 



now of fourth order. The new decay scales are given by 
A+^ = -i(^(^').n)|^(^')| 



A+(, = -i(fc(^) .n)|fc(^)| 
AT^ = -i(fc(^).n)|fc(^)| 
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(31) 



(32) 



(33) 



(34) 



with the abbreviations C = C" {q^) and S = S{qo). 
Hence we have for the real parts the relation 5R(A^^) = 

^(^tb) = -^(^j,a) = -^(^J,b) > 0' which means' that 
the solutions with superscript + are growing solutions 
and the ones with — are decaying solutions, all with the 
same range. Notice that they all also have oscillatory con- 
tributions, i.e. a non-vanishing imaginary part, ^(A) 7^ 0. 
We also have obviously = — A^* and A^^ = — Aj^. 
These relations are important for the proper matching of 
incoming and outgoing waves in the interface region be- 
tween the two grains. They imply ^(A^^) = ^{XJ^) and 
^(A^^) = ^(AJ^), therefore the oscillation frequency is 
equal for corresponding decaying and growing solutions. 
Also, the growth rates A^^ are only weakly imaginary in 
contrast to A^^. It is therefore not surprising that we find 
numerically that amplitudes of the strongly oscillatory 
solutions are very small, |c^^| <C \cf^\^ since an interface 
should mainly be a decay and not an oscillation of the 
density waves - the latter corresponds to a local change 
of the lattice spacing; the oscillatory modes can there- 
fore safely be neglected. For S iron we have e = 0.0860 
(i.e. e = 0.0923), and for this value A is only very slightly 
changed, and the amplitudes almost undistinguishable. 
Notice, however, that A formally becomes complex and 
that the decay rates 3?(A) for the incoming and outgoing 
waves are slightly diflFerent, but for present small values 
of e this diflFerence can be neglected. 



(i) - /X- A , - r\- \ IV. INTERFACE INTERACTION 

+ S> ^^P(\ta^) + exp(A+^2;) , (30) A. General framework 



with four independent solutions, since the equation is 



The simplest case of interacting solids is that of two 
lattices of the same material and with the same struc- 
ture that are perfectly aligned up to a translation in the 
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FIG. 6. Sketch of the geometry for shifted crystals. The 
displacement in the out-of-plane direction, Ay, is not shown. 
We assume that in the normal direction the crystals are not 
shifted, as illustrated by the dotted circles. 



contact plane, i.e. without misorientation between them 
(see Fig.[6|. If the crystals are fully aligned, which means 
that the atomic planes match, the interaction between 
the solid-melt interfaces is attractive, because exactly at 
the melting point (T = Tm) merging of the crystals re- 
moves two solid- melt interfaces, which reduces the total 
energy. On the other hand, the situation is more compli- 
cated if the crystals are shifted against each other, which 
implies elastic deformations of the lattices close to the 
grain interface. As we will show, a sufficiently large mis- 
match can lead to repulsive interactions. 

It is quite remarkable that the asymptotic of this struc- 
tural interaction between the crystals can be calculated 
fully analytically from the free-energy expression. The 
procedure is as follows: We assume two parallel crystal 
surfaces (see Fig.|6|, which are separated by a melt layer 
of thickness W. For large the density wave ampli- 
tudes are almost decayed in the center of the melt, and 
it is therefore sufficient to consider only the free energy 
terms up to quadratic order in the amplitudes. The re- 
lated equilibrium equations are therefore linear and can 
be solved easily, and the corresponding (approximative) 
solution has to be matched to the exact solution of the 
full problem of a localized interface at z ± W/2. By the 
means of this matching, we get an analytical expression 
for the disjoining potential. 

We start the analysis with the derivation of a conser- 
vation law. As before, we first ignore the higher order 
correction that stems from the box operator. The full 
free-energy expressions ([3| and ( 12 ) have the structure 



/ 



ifp + fk)dz, 



(35) 



where fp depends only on local terms (no gradients of the 
amplitudes), whereas fk contains only first order deriva- 
tive terms. Notice that due to the parallel structure, all 
amplitudes depend only on the coordinate z perpendicu- 



lar to the interfaces. Equilibrium demands 



(36) 



for all fields u^^\z). 

For a solid-melt-solid layer system, the free energy is 
in the spirit of equation ([T]) 



-WAf + V{W) + 27, 



(37) 



in the present case of the underlying NVT ensemble. The 
bulk free energy density difference Af = L{T — Tm)/Tm 
for a temperature deviation from the melting tempera- 
ture Tm corresponds to —AG introduced in Eq. ([T]) and 
will be discussed in detail below. 

To emphasize the analogy to a problem in classical 
mechanics, we use a dot for the spatial derivative in z 
direction. The "Hamiltonian" , 



H = fk - fp 



(38) 



is then a "constant of motion" , i.e. it does not depend on 
the z coordinate, 



H = 0. 



(39) 



For interfaces that are far apart, the amplitudes have 
almost decayed to zero in the melt region, and all contri- 
butions which are higher than qudratic in the amplitudes 
give only negligible corrections. Therefore, we get 



N/2 ^ 



Siqo) 

J 

^\c" {qQ)(k^^^ • nfu^^^u^^^'^y (40) 

The corresponding linearized "equations of motion" 
which describe the small amplitudes in the liquid region 
are therefore again given by Eqs. (23)-(25), with the only 



difference that we have here two interfaces, and therefore 
both exponentials are present. From this solution we can 
calculate the Hamiltonian in a quadratic approximation. 



H -2nokBT 



^ N/2 



(41) 



We can choose the origin z = in the center between 
the two interfaces, and then the exponential prefactors 
have the same absolute value but can differ by their 
phase, Cj^out = Cj^jn exp(z^^). Furthermore, from the 



general solution (24) it is obvious that a translation of 
the interface position in the normal direction increases 
or decreases the prefactors Cj by an exponential factor 
exp(AjAz), where Az is the shift distance. Therefore, we 
get Cj^in = Cj^o exp(— Aj W/2); the matching constants Cj^o 
were determined already in the previous section. Hence, 



1 



N/2 



\cjfi\'^ exp{-\jW) cos 4>j. (42) 
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In general, it is necessary to introduce a tilt term to fa- 
vor either the liquid or the solid state, because otherwise 
a repulsive or attractive interaction between the inter- 
faces would forbid the existence of a stationary solution 
(stable or unstable) with a specific melt layer thickness 
W. We therefore have to raise or lower the free en- 
ergy density of the solid phase relative to the liquid by 
A/ = L{T — Tm)/Tm' In particular, overheating above 
the melting point corresponds to A/ > 0. Notice that 
for the following calculation of the asymptotic interface 
interaction the precise form of the coupling in is not 
important, and only the tilt A/ enters into the result, 
provided that the bulk states u^^^ = and u^^^ = Us 
are temperature independent. This is the case e.g. for 
the coupling function (19), or — more generally — if 



which follows directly from Eqs. (42) and (43). This ex- 



the coupling function does not have a linear term in the 
amplitude variation 8u^^^ = u^^^ — Us in the solid and 
(5^/0) = 1^0) in the melt phase. The case that a linear 
term exists will be discussed in more detail below in Ap- 
pendix [Aj 

In the solid, the amplitudes are (up to a phase factor) 
all equal to 2X5, the gradients vanish, and therefore the 
Hamiltonian becomes 



H = -A/. 



(43) 



Comparison of this exact value, calculated from the solid 
phase, and the asymptotic value for large interface sepa- 
rations ( [42| ), calculated from the liquid phase, using the 
conservation law ([39| we obtain an implicit relation for 
the (asymptotic) width of the liquid layer W as function 
of the deviation from the melting point. A/. Asymp- 
totically, only the slowest decaying density waves with 
the smallest inverse decay length Xmin contribute to the 
Hamiltonian, and we get 



^,o\'^ ex.p{-XminW) COS ( 



-L 



T-Tm 



thus 



W 



In 



L{T-Tm) 



4no/CB^|Cmm,OpCOS( 



(44) 

This expression diverges logarithmically at the melting 
point, where ly = 00 is the equilibrium solution. If 
cos ^rnin IS positivc, wc find an asymptotic solution only 
for T > Tm' The interfaces attract each other, and this 
has to be compensated by overheating, i.e. favoring the 
liquid phase. On the other hand, for cos (j)min < 0, we 
have repulsive solutions asymptotically only below the 
melting point. 

At shorter distances the other density waves also con- 
tribute, and we therefore have to sum over all of them, 
which leads to an implicit relation for the melt layer 
thickness as function of temperature. 



N/2 



T-Ti 



M 



-M 



pression is valid as long as the overlap of the density 
waves is still small, such that the nonlinear energy con- 
tributions can be neglected. 

We can interpret the free energy shift A/ as the chem- 
ical force that balances the interface interaction. In fact, 
for a single interface it is the driving force for melting or 
solidification. From the equilibrium condition F'{W) = 
we get by means of Eq. ( [37| ) 

(46) 



H = -V\W) = -Af, 



where —V\W) is the disjoining force, which is derived 
from the disjoining potential ^(W). Integrating there- 
fore gives 



V{W)^-2nokBTJ^^S^x 



S{qo) 



N/2 



X "^\k^^^ ■ h\\cj ^of COS (I) jexp{-XjW), (47) 
3 

where we normalized the potential such that it decays to 
zero for infinitely far separated interfaces, in agreement 
with Eq. (37). Eq. (47) is the central result of this ar- 



ticle. Notice that the above expression of the disjoining 
potential is valid asymptotically for 1^ ^ 00. In this 
limit WXj ^ 1, which means that the interface thickness 
is small in comparison to the grain separation W. Then 
the interfaces are sharp, and the melt layer thickness is 
(uniquely) well-defined. For shorter distances, we use the 
same measure for W as defined above in Eq. (27), taking 



(45) 



into account that for the shifted crystals the interfaces 
remain planar (the amplitudes depend only on z). 

The solution of the linear equations is only valid for 
"non-slaved" fields, in particular those with A;^-^^ • n 7^ 0, 
and therefore these fields contribute differently to the 
disjoining potential by higher order nonlinear it ies. How- 
ever, in the above expression (47), the fast decaying and 
therefore negligible fields do not contribute due to or- 
thogonality, k^^"^ ■ n = 0. 

We can choose the origin of the coordinate system such 
that the amplitudes of one crystal are purely real. We as- 
sume that the other crystal is translated against it in the 
plane of the grain boundary, so (for a three-dimensional 
system) we have two translational degrees of freedom. 
The translation vector, Ar then obeys Ar • n = 0, so 
with h = z we get Ar = Axx -\- Ayy. The original non- 
shifted crystal is described by the expression 

5n{r) = no'^u^^\r)exp{ik^^^ • r), (48) 

3 

and a translation is therefore described by 

5n{r) = no ^ 'U*^-^^ (r) exp [ik^^^ -{r^Af)] 
3 

= no ^ u^^\f) exp[z^*^-^^ • Af] ex.p{ik^^'^ • r). 
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The complex shift factors are therefore given by 

(t)j = k^^^ • Ar. (49) 

We define the lateral dependence of the disjoining poten- 
tial for the fields with equal decay length A, or e quiv a- 
lently the same value /c • n, in agreement with Eq. (47) 



N/2 

cos^j-, 

j,k^^^ •h=k-h 



(50) 



where we sum over all amplitudes j with equal decay 
length. 

All density waves with the same decay length, i.e. equal 
value oi k-n and \j = A, have the same exponential decay, 
and we can define 



Vf^,^{W) = 2nokBT^^^^\k ■ h\ exp(-Al^). (51) 

The disjoining potential therefore becomes a superposi- 
tion of terms which factorize into a interface separation 
and translation part, 

N 

V{W, Ax, Ay) c^"YS{k- n, k^^^ • n)] x 
j 

\c,,o\'fU^x,^y)Vf^-niW), (52) 

B. {100} surfaces 

We consider two parallel (100) interfaces of bcc crys- 
tals as a first example. Asymptotically, the interac- 
tions stem exclusively from the density waves with the 
slowest decay, i.e. with the highest directional cosine 
k^^"^ ■ n. In this case, the principal reciprocal lattice vec- 
tors [110], [101], [110], [101] (and their inverses) have the 
same decay length, and the remaining, [Oil], [Oil] (+ in- 
verses) form a second group. All density waves within 
the same group have the same absolute value, but usu- 
ally differ in phase; notice that the amplitudes depend 
on the lattice shift. This is shown in Fig. [7[ where the 
absolute value of the density wave amplitudes is plotted 
as function of the position normal to the interfaces for 
a case without lattice shift. In the solid, all amplitudes 
reach the same bulk value Us due to the crystallographic 
symmetries. 

In agreement with the notation of (100) interfaces, we 
use a coordinate representation for n^^^^) = ^(1^0) = 
(1,0,0), and hence the tangential vectors have the co- 
ordinate representation x^^^^"^ = (0,1,0) and = 
(0, 0, 1). The translation is periodic with respect to shifts 
by one lattice unit a = 2V^7r/g'o in each direction x and 
y (the factor ^/2 comes from the fact that the recipro- 
cal lattice vectors point along the face diagonal of the 
bcc crystal). We can therefore introduce rescaled coor- 





-0.5 0.5 

[110], [101], [liO], [loi] — 
[oii],[oii] --- 
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FIG. 7. Absolute value of the density wave amplitudes for 



(100) 



0. 



(100) interfaces without lattice shift, ip^ 
The inset is a magnification around the origin, showing that 
the absolute value of the amplitudes varies smoothly there. 
The slaved fields decay quickly and do not show a visible 
overlap of incoming and outgoing waves. 





k<-^^ ■ n 


Phase shift 0^ 


Matching constant Cj,o 


[110] 


1/V2 


/ (100) 


0.165 


[101] 


1/72 




0.165 


[110] 


1/V2 


/ (100) 


0.165 


[101] 


1/V2 




0.165 


[Oil] 





(slaved) 


(slaved) 


[Oil] 





(slaved) 


(slaved) 



TABLE I. Matching and shift properties of the density wave 
amplitudes for (100) interfaces. 



properties are then 27r periodic for this coordinate repre- 
sentation. 

For the (100) interfaces, the phase shifts of the fields 
are summarized in Table [T| together with the previously 
determined matching constants, and we therefore obtain 
for the longest range exponentials 



.(100) 



(100) 



,4^»)) = -2(cos^r^ 



• COS '02"'^^^^ 



), 

(53) 

which is plotted in Fig.[8| The disjoining potential for the 
most attractive situation (matching lattices, (a)) and the 
most repulsive case (b) are shown in Fig. |9] The predic- 
tions are compared to numerical calculations, which were 
obtained in a dynamical run at T = Tm- This means that 
we set up a solid-liquid-solid "sandwich" structure, with 
a phase shift between the solid phases in the real space 
implementation. Due to the overlap of the interface pro- 
files we have an attractive or repulsive interaction be- 
tween the interfaces, thus the configuration is not in full 
equilibrium. During the time evolution we numerically 
compute the melt layer thickness W and the energy F 
according to Eq. (12), without the correction term from 



dinates = Axqo/\^ and ip. 



(100) 



Ayqo/V2; ah 



the box operator. The dependence F{W) is then plotted 
and compared to the analytical predictions. This method 
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y/^°°' = [010] 



W'°°' = [001] 



FIG. 8. Lateral dependence fi-\i^ slowest decaying 

contribution to the disjoining potential of two parallel (100) 
surfaces. The interaction is repulsive for positive function 
value. The two arrows mark the most attractive (a) and most 
repulsive situation (b). 



> 




7 



10 



FIG. 9. The disjoining potential for the two cases (a) (solid 
curve) and (b) (dashed curve) of (100) interfaces. For each 
case, the squares show the result from the numerical simula- 
tion, the lines the asymptotic prediction, taking into account 
the slowest decaying density waves. 



is approximative in the sense that the system is not in 
full equilibrium with dA^^^ /dt = 0. A more precise ap- 
proach is to balance the interaction with a thermal tilt 
T 7^ Tm and then to calculate the energy for the relaxed 
solution F — Ft] this approach is used for the interaction 
of misoriented crystals, which are treated in Section [V[ 
As long as the interaction is weak, both methods give the 
same results, and we have checked that the present results 
are robust. Also, they agree very well with the analyt- 
ical predictions for the asymptotic interaction. For the 
special case that the disjoining potential has a minimum 
- a case that we will encounter later -, the dynamical 
runs converge to this point, where the interaction energy 
therefore becomes exact. 

In the liquid region, the density wave amplitudes are 




3.5 

3 [ 
2.5 

2 
1.5 

1 

0.5 


-0.5 



[110], [101], [110], [101] 
[Oil], [Oil] 



-10 
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FIG. 10. Top: Absolute value of the density wave amplitudes 
for (100) interfaces for crystals that are shifted by half a lattice 
unit, ij:^^^^ — i^^^^^ — 7T. The inset is a magnification around 
the origin, showing that the absolute value of the amplitudes 
has a cusp there. The slaved fields decay quickly and do 
not show a visible overlap of incoming and outgoing waves. 
Bottom: Corresponding phases of the amplitudes. 



given by the expression ( 24 ) . It is instructive to look also 
at amplitude and phase separately. With a real coefR- 
Cj,out = Cj^in exp(z0j) we obtain at z = 



|zi(^-)|2 = 24. „(1+ cos 



(54) 



For the special case = tt (the most repulsive case) 
the amplitudes have a cusp there, and correspondingly 
the phase jumps discontinuously. Notice, however, that 
this singular behavior appears only in the polar repre- 
sentation of the complex amplitudes; in a complex sense 
they are smooth at 2; = 0. This behavior is visualized in 
Fig. [To) We note that a cusp in the "order parameter" 



was introduced phenomenologically in Ref. [24J Here it is 
a natural consequence of the description. 



C. {110} surfaces 

The situation immediately becomes more complex for 
the next example of (110) interfaces, thus n^^^^) = 
(1,1,0) /a/2. The tangential vectors are here defined 
through the coordinate representation x'^^^^^ = (0, 0, 1) 
and = (1, — 1, 0)/a/2, and the phase factors are 
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v/''°' = [001] 




FIG. 11. Lateral dependence fi'^^^\ ,^ of the second contribu- 

tion to the disjoining potential of two parallel (110) surfaces. 
The interaction is repulsive for positive function value. The 
two arrows mark the most attractive (a) and most repulsive 
situation (b). 



^(110) _ Ax go /a/2 and il^2^^^^ = Ayqo/2^ to recover the 
27r periodicity. This means that the axes in the interface 
plane are stretched differently, and therefore the geome- 
try loses its fourfold symmetry (it reduces to a C2 sym- 
metry), in agreement with the fact that in the interface 
plane the distances between the atoms are different for 
the two perpendicular directions. 

Here, the density wave which corresponds to the prin- 
cipal reciprocal lattice vector k'^^^ = [110] has the longest 
range, since its crystallographic ordering extends the far- 
thest into the melt, k^^"^ ■ n = 1. Therefore, the longest 
range interaction is mediated by this density wave. Since 
it is a plane wave, it has no lateral dependence, which 
means that the disjoining potential does not depend on 
ip2^^^\ and this contribution to the interaction is 
always attractive. This implies that at large distances we 
always find an attraction of crystals, irrespective of the 
lattice shift, i.e. 



(110) 
k-h—1 



-1. 



(55) 



As soon as the interfaces come closer to each other, 
the contribution from the next density waves becomes 
noticeable. In this case, it comes from the density waves 
with reciprocal vectors [101], [Oil], [101], [Oil] (and 
their inverse vectors), which have all a directional co- 
sine k - n = 1/2, thus the range of their contribution 
to the interaction is only half of the range of the lead- 
ing term. The total lateral dependence from this set of 
density waves. 



(56) 



is shown in Fig. [TT] It has attractive and repulsive re- 
gions: If the crystals are perfectly aligned, V^^^^^^ = 





k^'^ ■ h 


Phase shift (pj 


Matching constant Cj,o 


[110] 


1 





0.116 


[101] 


1/2 




0.372 


[Oil] 


1/2 




0.372 


[101] 


1/2 




0.372 


[Oil] 


1/2 




0.372 


[110] 





(slaved) 


(slaved) 



TABLE II. Matching and shift properties of the density wave 
amplitudes for (110) interfaces. 
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FIG. 12. Disjoining potential for (110) interfaces. At large 
distances, the interaction is always attractive and solely deter- 
mined by the (110) density wave. The graphs shows numeri- 
cal results for two different shifts together with the analytical 
predictions: The dotted curve takes into account only the 
longest-range exponential, the solid (a) and dashed curve (b) 
also the corrections due to faster decaying density waves. 



ip2~^^^ = 0, the interaction is of course attractive, be- 
cause the interfacial energy would vanish completely if 
the crystals merge. For maximum mismatch, i.e. if the 
crystals are shifted by half a lattice unit in one direction, 
the interaction has reached the strongest repulsive situa- 
tion. For a shift by half a lattice unit in both directions 
we recover again the attractive case, because then the 
crystallographic planes in the (110) surface match again. 

The relevant data for the calculation of the asymptotic 
disjoining potential is summarized in Table [llj and the 
potential is plotted in Fig. [12] for the case of no misfit 



(a) and the most repulsive case (b) with ipi 



(110) 



TT and 

0, as illustrated by the arrows in Fig.|ll| As 
mentioned before, we have a purely attractive behavior at 
large distances independent of the lattice translation, and 
the agreement with the analytical prediction is confirmed 
in the logarithmic plot Fig. [13) At shorter distances, the 
numerically calculated disjoining potential deviates from 
the analytical prediction from the slowest decaying den- 
sity waves only (dotted lines in Fig. 12), and the inclusion 
of the next terms leads to a significantly better agree- 
ment (solid and dashed line) , and we observe the distinc- 
tion between the attractive and repulsive cases. Only 
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FIG. 13. Asymptotics of the disjoining potential for (110) 
interfaces. At large distances, the interaction is always at- 
tractive and solely determined by the (110) density wave. 
The graphs shows numerical results for two different shifts 
together with the analytical prediction (only the contribution 
from the slowest decaying exponential). 



at short distances is the interaction strongly affected by 
nonlinear contributions. In particular, we observe a sta- 
ble minimum in the disjoining potential for the most re- 
pulsive case, because a hard core repulsion due to elastic 
deformations prevents a full merging of the interfaces. 



D. {310} surfaces 







n 


Phase shift 0^ 


Matching constant Cj,o 


[110] 




i 0.89 




0.128 


[101] 


3/2V2-r 


^ 0.67 


^(310)^^(310) 


0.18 


[101] 


3/2^2 ? 


^ 0.67 


^(310) _ ^(310) 


0.18 


[110] 




i 0.44 




0.425 


[Oil] 


1/2^ 


^ 0.22 


(slaved by cubic) 


(slaved by cubic) 


[Oil] 


1/2^;: 


^ 0.22 


(slaved by cubic) 


(slaved by cubic) 



TABLE III. Phase shifts and matching constants for the (310) 
interface normal for the non-slaved fields. Notice that the 
fields with the highest scalar product hS^^ • n have the longest 
range. 




As a last example we investigate (310) interfaces, 
where we receive nontrivial contributions from the first 
and the second exponentials. With fi^^^^^ = (3, 1, 0)/ a/TO 
we use tangential vectors x'^^^^^ = (0, 0, 1) and y^^'^^^ = 
(1, —3, 0)/aAO- The 27r periodic in-plane coordinates are 
^(310) _ ^ ^ V^2^"^^^ = Ay qo / >/20, and all data 
is summarized in Table |III[ The longest range density 
wave is i^no, and the lateral dependence of the disjoining 
potential 



(310) 



- COS 2?/;2^"^^^ 



(57) 



is shown in Fig. 14 Two fields, [101] and [101] (plus in- 
verse vectors), contribute to the next exponential, there- 
fore the lateral dependence of this term is 



f - - cosf^^^'^^ 



- V^f^^^) -cos(V^i' 



(310) 



-2 cos ip^^^"^ cos ip^^^"^ , 



(58) 



see Fig. 15 



We investigate in particular three different shifts, all 
with ^ ^ 4^^°)'" = 0.04 • 27r, 4^^°)'^ = 

0.07 •27r and 7/^2^^^^'^ = 0.11 •27r. For this inclination, both 
the longest range term and the next term have - depend- 
ing on the mismatch - attractive and repulsive regions. 



The three scenarios are indicated by the arrows in Fig. 14 



FIG. 14. Lateral dependence of the longest range contribu- 
tion of the disjoining potential for the (310) interfaces. The 
potential is attractive in regions where the value is negative. 
In particular, this contribution to the interaction does not 
depend on the translation in the [001] direction. 




71/2 



¥2 



(310) _ 



[130] 



FIG. 15. Lateral dependence of the second longest range con- 
tribution of the disjoining potential for the (310) interface 
normal direction. The potential is attractive in regions where 
the value is negative. 
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FIG. 16. The disjoining potential for (310) interfaces for the 
three cases a, b, c, as explained in the text. For each case, 
the isolated points show the result from the numerical simu- 
lation, the dotted lines the asymptotic prediction, taking into 
account only the slowest decaying density wave [110], and the 
solid line the analytical predictions, using the first and second 
exponentials. The potential is here in all cases attractive at 
large distances and repulsive for small separations W. 



Analytical asymptotics 



Numerics □ 



> 




10 11 12 13 



FIG. 17. Logarithmic plot of the asymptotics of the disjoining 
potential for V^f = tt and V^^^^°^'"^ = 0.04 • 27v for (310) in- 
terfaces. The numerical results are compared against the an- 
alytically determined asymptotic behavior including the two 
slowest decaying exponentials. 



for the lateral dependence of the slowest decaying expo- 



nential and also in Fig. 15 for the next exponential. Ap- 
parently, for the sample cases (a), (b), and (c) the slow- 
est decaying exponential is always attractive, whereas 
the second is repulsive. From (a) to (c) the strength of 
the first exponential becomes smaller, and therefore we 
see a crossover from a long-range attraction to an inter- 
mediate repulsion. This prediction is confirmed by the 
numerical results in Fig. [iGj Fig. [Tt] shows the asymp- 
totics of the disjoining potential for = tt and 
^(310), a _ Q 2^ ^j^^ ^i^Q comparison with the analyt- 
ical prediction. 



E. Box operator corrections 

Let us briefly discuss the influence of the correction 
terms of the box operator, which have been neglected 
in the discussion so far. With the previous knowledge 
that the additional terms are small, it is immediately 
transparent that the results can be modified only slightly, 
and therefore the above simplified picture remains valid. 
Nevertheless, the analysis can also be formally performed 
in this more complicated case, and this is outlined here. 

The presence of the box operator leads to the following 
modifications: First, the amplitudes in the linearized re- 
gion become a superposition of four exponential solutions 
instead of only two. Determining the disjoining potential 
now requires identifying matching pairs of incoming and 
outgoing waves in the sense of a conservation law. Two 
of the exponentials are strongly suppressed, since they 
show relatively fast oscillations. Second, the concept of 
the Hamiltonian as in classical mechanics is only appli- 
cable if the free energy density contains only first order 
derivatives (the kinetic term). The box operator, how- 
ever, introduces higher order derivatives, and therefore 
this concept has to be generalized. 

A generalized conservation law, which is valid also for 
misoriented grain boundaries, where the box operator is 
essential, is derived in Appendix [Bj and here we need 
only the special case that all fields depend only on the 
coordinate normal to the grain boundary. (In the gen- 
eral case, the Hamiltonian is an integral expression along 
the grain boundary plane, which reflects the fact that 
the interaction forces can vary spatially and have to be 
averaged to get the net force.) Here, the interaction is 
homogeneous, and therefore the following expression be- 
comes a conserved quantity: 

N/2 



(59) 



where / is the free-energy density. It corresponds to a 
generalized Legendre transformation, with "momenta" 



P 



df 



df 

dm ' 



(60) 



where we treat u''^^ and the complex conjugate i/^-^^* as 
independent functions. The Hamiltonian is conserved, 
i.e. H = 0. From the linearized solution (30) we obtain 



after some straightforward but tedious algebraic manip- 
ulations 

TT _ "^n^kBT ^ / _^ + _ _^ + _ 

(61) 

for the value of the Hamiltonian, calculated in the liquid 
up to second order. 
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Similar to before, the prefactors of the exponentials in 
the linearized solution decay with the melt layer thick- 
ness, and we have 



= c-;%xp(A-,iy/2), 



clfexp{-X+W/2), 



(62) 
(63) 
(64) 
(65) 



Therefore, the disjoining potential becomes 



V{W) 



ZnokB-L 



j 



-,0* +,0 

c c 



exp(-A+ ly) 



,0* +,0 



exp(-A+,l^) 



c.c. 



(66) 



Again, the prefactors acquire a complex factor if the crys- 
tals are translated against each other. As mentioned be- 
fore, the prefactors cjj^ and c^f are small and can be 
neglected. For rough interfaces (e 0) we recover the 
above expression (47) for the disjoining potential. 



V. INTERACTION BETWEEN MISORIENTED 
GRAINS 

An analysis for the long range interaction as for the 
shifted crystals is not possible if a misorientation is in- 
volved, since the problem is not one-dimensional any- 
more. However, the range of the interactions can still 
be understood using similar arguments, and the central 
outcome is that they are significantly shorter ranged. We 
consider here the case of a tilt grain boundary, to illus- 
trate the basic idea. As had been discussed in detail in 
Ref.|40]the presence of a lattice rotation makes the use of 
the full box operator mandatory, and still the description 
is only valid for small misorientations. 

Let us assume that the left grain {z — oo, charac- 
terized by subscript — ) is rotated by <!>_ with respect to 
the reference orientation of the RLVs, whereas the right 
grain ^ oo, subscript +) is rotated by Still the 
interfaces are assumed to be planar and normal to the z 
axis. 

First, the left grain has amplitudes 



( exp [ifc(i)t]V[(<E>_)f] 



(67) 



in the bulk, where the dagger denotes transposition and 
M($) = R(^) — I with the identity matrix I and the 
rotation matrix R(^), 



R($) 



COS <l> sin $ 
— sin ^ cos ^ 



Similarly, for the right grain 



(68) 



(69) 



This suggests looking for solutions in the liquid region of 
the structure 



c-^ exp[i^(^)tM($_)f]exp (A-J$_)z) 
cj^^ exp[ik^^^'^M{^_)r\ exp (^Xj^^{^_)z^ 
exp[i^(^-)tM($+)r] exp (A+,($+)^) 
c+,exp[i^(^)tM(^+)r]exp (a+,($+)z) , (70) 



in analogy to Eq. (|30|), with ^{X~J < 0, 3?(A^"J < 0, 



3?(A+^) > 0, 3?(A_^5) > 0. The ranges ^f^/i, are computed 
from the linearized equilibrium condition (28) with the 
help of the rotation theorem 



□I 



/(f) exp(ip)1'Mr)l = exp(ifc(^')^'Mf)n2+/(f), 



for any function /(r) and 



2qo 



(71) 



with the rotated reciprocal vectors = R'*'^'^-^^ see 
Ref. |40] for details. It turns out that the decay ranges 
are given by the same expressions as for the shifted crys- 
tals, Eqs. (31)-(34), but the reciprocal vectors have to be 



rotated here appropriately in the k^^"^ ■ n term. 

Inserting these expressions into the generalized conser- 
vation law (B6) using only quadratic terms gives H = 
in disagreement with the tilt H = —A/. This shows 
that the longest range interaction is not mediated by 
the quadratic terms in the functional but stems from the 
higher order nonlinearities. Since their contribution van- 
ished quickly in the melt phase, it is intuitively clear, that 
the interaction range for misoriented grains is shorter 
than for shifted crystals. 

One can interpret this statement also in a physical 
way: For two misoriented grains the normal shift between 
lattice planes of the two crystals varies along the grain 
boundary, and therefore the interface alternatingly con- 
sists of region, where the atomic planes match and where 
they are out of phase. This leads to alternations of at- 
tractive and repulsive regions along the grain boundary. 
Since the strength and size of attractive and repulsive 
regions is equal at quadratic order, their contributions 
cancel each other in the total interaction energy between 
the grains. Thus, only shorter range higher order terms 
can be responsible for the disjoining potential here. 

To understand the range of the remaining interaction 
further, one can continue to employ the conservation law 
(B6). The next step is to assume that the interaction 



stems from the cubic terms in the functional. They can 
appear in two different ways: First, products of two in- 
coming and one outgoing waves or one incoming and two 
outgoing waves from expression (70). Second, the cu- 



bic nonlinearities (which appear as quadratic terms in 
the equilibrium conditions) generate perturbations of the 
basic solution ( [TQ] ). The structure of these perturba- 
tions du'^^^ would be again be a product of two density 
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FIG. 18. Comparison of the interaction decay for (100) shifted 
crystals and a (100) symmetric tilt with 20 — 45.2° for ^-iron. 



waves, and a product of the type u^^^5u^^^ would then 
have then same structure of a product of three density 
waves. However, such a product would contain a terms 
like exp [(Ai + A2 + As)^] with three decay lengths Ai, A2, 
A3, which would appear in the conservation law. Since it 
has a nontrivial z dependence, all these terms in the end 
have to cancel, since by the conservation law the Hamil- 
tonian must be z independent. Therefore, also the cubic 
terms cannot contribute to the long-range interaction. 

At quartic order, terms like exp [(Ai + A2 + A3 + \4)z\ 
can again be z independent in the Hamiltonian, and they 
do not cancel, thus they finally balance the temperature 
term A/. The interaction range is therefore set by the 
sum of two A values, and decays therefore twice faster as 
for shifted crystals. This is also confirmed by numerical 



investigations. Fig. 18 shows the cases of the shifted crys- 
tals in comparison to a repulsive grain boundary. The 
shifted crystals case depicts the maximum repulsion for 
(100) interfaces, as discussed before in Fig. [9] (curve b), 
again without the correction term of the box operator. 
In comparison, the disjoining potential for the repulsive 
grain boundary with a misorientation of 20 = 45.2° de- 
cays significantly faster. For these simulations we used 
e = 0.0923 and also took into account the box operator 
correction, since otherwise rotated crystals would melt 
spuriously at T = Tm- For the mapping to physical data 
we used the parameters for bcc a-iron, which were pre- 
viously determined in Ref. i38| 



Fig. [19] shows the same data on a logarithmic scale, 
together with the analytical prediction of the slowest de- 
caying quartic interaction term, which stems from the 
[110] density waves; here we note that for the prediction 
of the decay range also the rotation of the interface nor- 
mal by has to be taken into account in expression (25). 
The prefactor of the interaction energy is matched to the 
computed disjoining potential. These numerical results 
confirm that the interaction of misoriented crystals is me- 
diated by quartic terms in the framework of this model, 
and therefore the interactions of wet grain boundaries is 
very short ranged. 




0.01 



0.0001 > 



1e-06 



FIG. 19. Comparison of the interaction decay for (100) shifted 
crystals and a (100) symmetric tilt with 20 = 45.2° for ^- 
iron. The decay range of the quartic asymptotics is calculated 
without including the box operator correction term, which is 
negligible here. 



VI. DISCUSSION AND SUMMARY 

We have calculated the interaction between solid-liquid 
interfaces based on amplitude equations, which are de- 
rived from PFC or density functional theory. In the 
framework of this model, the tail of the structural interac- 
tion can be calculated fully analytically for grains which 
are not misoriented but only differ by a lateral transla- 
tion. It is short-ranged and decays exponentially with 
the grain separation. Depending on the lattice mismatch 
we find that the interaction is either attractive or repul- 
sive. It is most attractive if the lattice planes are fully 
aligned, such that complete freezing would remove any 
interface between the crystals. In the opposite extreme 
case, that the grains are shifted against each other, such 
that a strong mismatch appears when the liquid layer dis- 
appears, leading to strong elastic deformations, the inter- 
action is repulsive. The entire interaction is a superpo- 
sition of the contributions of the different density waves, 
and the longest range fields at a solid-liquid interface, 
i.e. those density waves which extend the crystalline or- 
dering furthest into the melt, also give the longest-range 
contribution to the solid-melt interface interaction. The 



range of this interface interaction is given by Eq. ( 25 ) for 
the individual density waves. This analytical expression 
also shows that the range of the interaction is determined 
by scattering properties of the melt phase and the rela- 
tive orientation of the density wave vector to the interface 
normal. 

The case of purely translated grains describes unsta- 
ble configurations or repulsive interactions, since strong 
lateral forces will force the system back to configura- 
tions with aligned crystallographic planes. Nevertheless, 
the results shall be relevant for the understanding of 
the merging of dendrite sidearms from the same grain, 
where due to elastic deformations the lattices in the side- 
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branches are shifted against each other. Apart from the 
remarkable feature that the interaction and its origin can 
be understood fuhy analyticahy, we also mention the re- 
lation to the similar concept of 7-surface&^^ - the energy 
landscape of two half crystals are which are tangentially 
displaced against each other - which is essential for the 
understanding of stacking faults and other defect forma- 
tion mechanisms. Here, we obtain a fully analytical pre- 
diction of this lateral dependence of this energy landscape 
for solid- melt-solid layer systems at larger distances. 

At this point we also mention that the term interac- 
tion must not be interpreted as a mechanical force that 
leads to an interaction between the grains, but rather 
a thermodynamic force. The central difference is that a 
mechanical attraction e.g. would move the grains towards 
each other, so each atom moves. Here, however, we con- 
sider the situation of a melt layer that separates the two 
grains, and therefore the attraction between the grains 
would manifest in the solidification of the melt layer. As 
a consequence, the gap between the crystal is closed, but 
without a rigid body motion of the entire grain. In other 
words, during the solidification process the number of 
atoms in the solid phases increases, whereas it would be 
conserved for a purely mechanical motion. This aspect 
is also important from another point of view: In the con- 
sideration of the shifted crystals we have excluded trans- 
lations in normal direction {z direction) and only inves- 
tigated motion in the tangential xy plane. This means 
that the atomic planes are always aligned in the normal 
direction and only exhibit a mismatch in the others. One 
could also consider the translation in z direction, which 
would then lead to an additional oscillatory interaction 
dependence in this direction. This Az dependence is not 
related to the exponential decay of the disjoining poten- 
tial, which appears separately on a larger scale. 

Beyond the case of pure grain translation we have also 
considered grain boundaries with a misorientation. In 
this case, a full analytical calculation of the disjoining po- 
tential is not possible anymore and numerical investiga- 
tions are needec^. Nevertheless, we have explained that 
the interaction then stems from higher order terms in the 
free energy functional, since the longest-range contribu- 
tions from quadratic terms cancel. As a result, we find 
that the disjoining potential decays twice faster then for 
shifted crystals. This prediction of the interaction range 
is confirmed also by numerical simulations and further 
discussed in Ref. [311 



Appendix A: Linear temperature coupling 
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The difference between the coupling functions (18) and 



(20) is that the first is quartic in the amplitude varia- 
tions in the bulk states and the latter linear. From this 
we immediately conclude that in the first case it is not 
necessary to take the tilt term into account for the so- 
lution of the linearized equations in the liquid region, as 
it is of higher order, and it appears there only effectively 
through the energy shift as the value of the Hamiltonian. 
The situation is different for the linear coupling, and dif- 
ferent effects have to be taken into account formally: (i) 
The amplitudes in the liquid change since the minimum 
of the potential energy is shifted away from u^^^ = for 
T ^ Tm, (ii) the amplitudes in the solid change since the 
minimum of the potential energy is shifted away from 
u^o) = for T 7^ Tm, (hi) the thermal tilt gives a con- 
tribution to the Hamiltonian in the liquid, which needs 
to be taken into account up to second order, and (iv) the 
value of the Hamiltonian changes in the solid due to the 
shift of the solid amplitudes. 

At a first glance, one may therefore expect that the 
results are changed by these effects, and that the long 
range interaction should depend on the precise choice 
of the coupling function. However, here we show that 
this is not the case up to linear order in (T — Tm)/Tm- 
For simplicity, we again do not take into account the 
correction term from the box operator. 

Writing /({^x(^)},T) = fnoni{{u^'^}) + It in the spirit 
of Eqs. (|3| and ( [21] ), where fnoni therefore contains the 
free energy density terms to cubic and higher order (lead- 
ing to the nonlinear terms in the equilibrium conditions) , 
we therefore obtain the equilibrium conditions 



8F 



a/„o„/({w(^)}) , 5/t({w(^)}) 



0. 



(Al) 



The equilibrium conditions for the hquid between two 
sohd phases are therefore up to first order in the ampH- 
tudes 



noksT 



U(9o) 



+ L 



T-Tm 1 



'-M 



Nus \uii) 



= 



(A2) 



and differ from the previous condition (23) only by the 



temperature term. For simplicity, we consider from now 
only the case of real amplitudes, i.e. no crystal transla- 
tion, u^^^ /\u^^^\ = 1. Then the general solution is 

^(j) = c ^ Cj^in exp(-A_^-2;) + Cj^out exp(A_^-z) (A3) 

where the decay parameters are unchanged and given by 
Eq. (25). The constant term is 



S{qo) ^T-Tm 1 



noksT 



'-M 



(A4) 
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In the same way, we can analyze the behavior in the 
solid phases, where the amplitudes are constant. One 
readily finds that the deviation from the previous bulk 
value u^^"^ = Us is linear in the temperature deviation 
{T-Tm)/Tm. 

The expression for the Hamiltonian up to second order 
in the amplitudes also now contains the tilt term. 



H 



noksT 



(A5) 



Inserting the above solution gives up to first order in 
(T — Tm)/Tm the same expression ( [41] ) as before. 

In the solid, the value deviates from the previous ex- 
pression H = — A/ = L(T — Tm)/Tm only by quadratic 
corrections in A/. Consequently, the long range interac- 
tion is unaffected by the choice of the linear instead of a 
higher order coupling function for the temperature. 



Appendix B: The generalized conservation law 

As has been demonstrated in Section IIVI the consid- 
eration of the Hamiltonian i^" as a conserved quantity 
is a valuable way to understand the interaction of two 
crystals, which have the same lattice orientation but are 
shifted against each other. The limitations were (i) the 
neglect of the box operator correction and (ii) the re- 
striction to pure translations, which forbids the analysis 
of grain boundaries, where the crystals are misoriented. 

Here we generalize this concept to overcome these re- 
strictions. 

For the shifted crystals, the translation leads to a 
multiplication of the amplitudes by a spatially constant 
phase factor. The main difference is that a rotation leads 
to non-constant phase factors. In particular, the ampli- 
tudes then do not depend only on a single coordinate 
normal to the grain boundary, but on all spatial coordi- 
nates. Therefore, a proper conservation law should also 
take into account the directions parallel to the interface. 

We consider a "Lagrangian" (i.e. the free energy den- 
sity in the present context) of the form 

C = C{{u^^^}, {ii^'^}, {i^(^>}, {u^'^'}, {^x^^'^*'}, 

{ii^^)}, {i^(^>}, {u^'^"}, {u^'>"}) (Bl) 

where expressions like {u^^^} denote the set of all am- 
plitudes u^^\ For simplicity, we assume that all fields 
depend only on two coordinates, which we choose to be 
normal and tangential to the interface. This is the case 
e.g. for tilt grain boundaries, and a twist would require 
the straightforward inclusion of another tangential de- 
pendence. Derivatives with respect to these directions 
are denoted by a dot for the normal and a prime for the 
tangential direction, although this assignment of direc- 
tions is in principle arbitrary; however, it will turn out 



to be a useful choice. We align our coordinate system 
such that X is the normal and y the tangential direction. 
In contrast to the pure Hamiltonian system in the pre- 
vious section here also higher order derivatives appear. 
It turns out that in our case mixed mode derivatives like 
ii^^^ do not appear, and therefore we do not consider 
them. 

We introduce generalized momenta. 



dC 
dC 

dm' 



dC 
dC 



(B2) 
(B3) 



Here, u^^"^ and u^^"^* are treated as independent variables. 
The equilibrium conditions for the amplitudes 



jCdr = (B4) 



can be written as 
dC 



jUY + rU) + ^Or' = 0. (B5) 



We obtain then the following conserved quantity, which 
means H = 0: 



H 



N/2 



) 



c 



(B6) 



Here, N/2 is the number of independent density waves 
(we write the complex conjugate fields explicitly and 
must not double-count them). The proof is straightfor- 
ward: Performing the normal derivative and application 
of the equilibrium conditions (B5) yields after a few al- 
gebraic simplifications 



H= dyY^ -dy (^q^^^m + q^^>u^^>^ 

-dy (^S^^^m' ^S^^>U^^>'^ ^dy (^S^^y^^ 

= 0, 



where the last step follows from periodicity along the 
grain boundary (in y direction). 

The expressions are written here for a two-dimensional 
dependence of the fields, which is the case of tilt bound- 
aries in a three-dimensional system. It is obvious that 
for more general cases, e.g. twists, the concept can easily 
be generalized by introduction of a second coordinate in 
the grain boundary plane. 
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